1.1 The beginning

At a yelp-highly-rated chicken place one summer, I was waiting for my order, which took forever (exaggerated). Without any other better things to do, I re-studied the menu there, trying to understand the pricing structure. See photo:

I was wondering to myself tons of questions: how can I get the best deal? What is the cost per wing, and cost per drum? Does the restaurant use a formula to determine the pricing?

Bugger! “I am cursed. Shouldn’t have come here. I might not get out of here alive…” I thought at the time.

You can see the set prices in the photo. We will need to re-format our data, from this wide format (often in pivot and summary tables), to a flat, long format here. Let us practice our Markdown skill at the same time to tabulate them:

Let us now enter these as a dataframe in R.

# korean-fried-chicken-wing pricing
kfcw <- data.frame(wing=c(-5,10,15,20,40,0,0,0,0,3,7,12,20), drum=c(0,0,0,0,0,3,5,10,15,2,4,6,9), price=c(6.69,12.99,17.99,23.99,45.99,6.69,10.99,19.99,29.99,7.99,16.99,27.99,43.99) )
xkabledply(kfcw) 

So there are 13 data points in our dataset. (Notice the use of inline R codes here.)

1.2 Exploratory Data Analysis (EDA)

There are a bit of things we typically look at for EDA.

1. Basic statistics
   • mean, s.d., median, range (four spaces at the start of line for proper sub-list indentation)
2. Simple correlations and tests
   • correlation matrix if applicable
   • z-test, t-test, anova test if applicable
   • chi-squared test if applicable
3. Normality
   • QQ-plot
   • boxplot
   • histogram
   • Shapiro-Wilk test
   • …

# str(kfcw)
# summary(kfcw)
xkablesummary(kfcw) # better display than the line above

kfcw[1,1]=5  # equal sign "=" and assignment operator "<-" are interchangeable in R
# structure of the data frame kfcw
# str(kfcw)
# summary(kfcw)
xkablesummary(kfcw) # better display than the line above

loadPkg("corrplot")
corrmatrix = cor(kfcw)  # more detailed pair-wise correlation test can be obtained from cor.test(kfcw$wing,kfcw$price) etc
corrplot.mixed(corrmatrix, 
               title="Correlation Matrix for Chicken Meal Price",
               mar=c(0,0,1,0) # fixes the position of title
               )

qqnorm(kfcw$price, main = "Price Q-Q Plot", ylab="Price Quantiles ($)") 
qqline(kfcw$price)

qqnorm(kfcw$wing, main = "Wing Q-Q Plot", ylab="Wing-count Quantiles") 
qqline(kfcw$wing)

qqnorm(kfcw$drum, main = "Drum Q-Q Plot", ylab="Drum-count Quantiles") 
qqline(kfcw$drum)

# We will learn to use the more powerful ggplot soon, instead of this generic boxplot function
boxplot(kfcw, col=c("red","blue","green"), ylab="count or price($)", main="Boxplots for the three variables")
axis(side = 4)

# We will learn to use the more powerful ggplot soon, instead of this generic hist function for histograms
barcolors = c("green", "violet", "orange", "blue", "pink", "red", "yellow", "cyan")
hist(kfcw$price, main = "Histogram for Price distribution", xlab="Price ($)", col=barcolors, breaks = 10)

hist(kfcw$wing, main = "Histogram for Wing-count", xlab="Wing Count", col=barcolors, breaks = 6)

hist(kfcw$drum, main = "Histogram for Drum-count", xlab="Drum Count", col=barcolors, breaks = 6)

priceshapiro = shapiro.test(kfcw$price)
wingshapiro = shapiro.test(kfcw$wing)
drumshapiro = shapiro.test(kfcw$drum)

1.3 Linear Model

Now if those passed the smell test, we can try build some models. Linear model is first to come to mind. We will model the price with wing and drum as the two independent variables (also called features).

# build a simple linear model (least square fit) of price as a function of everything else.
chicklm = lm(price ~ ., data=kfcw)
# summary(chicklm)
xkabledply(chicklm, title = "Summary of linear model for Chicken Meals")

As we can see, the R\(^2\) (using \(\LaTeX\) formatting here) value of the linear model is 0.9957, which shows the prices are set with rather strict per-piece-structure.

If we use the (default 95%) confidence intervals of the coefficients as shown here:

coeffconfint = confint.lm(chicklm) 
# coeffconfint
xkabledply(coeffconfint, title = "Coefficient Confidence Intervals (CI) at 95%") # display better than the line above

loadPkg("faraway") # faraway library is one of them has a vif function
# VIF check the collinearity issues between different variables/features in a (linear) model
# vif(chicklm)
xkablevif(chicklm, title="Chicken Model VIFs")  # better display than the line above

# unloadPkg("faraway")

loadPkg("plot3D")
# reference from http://www.sthda.com/english/wiki/impressive-package-for-3d-and-4d-graph-r-software-and-data-visualization
# x, y, z variables
x <- kfcw$wing
y <- kfcw$drum
z <- kfcw$price
# Compute the linear regression (z = ax + by + d)
# chicklm <- lm(z ~ x + y)
# predict values on regular xy grid
grid.lines = 20
x.pred <- seq(min(x), max(x), length.out = grid.lines)
y.pred <- seq(min(y), max(y), length.out = grid.lines)
xy <- expand.grid( wing = x.pred, drum =y.pred)
z.pred <- matrix(predict.lm(chicklm, newdata = xy), nrow = grid.lines, ncol = grid.lines)
# fitted points for droplines to surface
fitpoints <- predict.lm(chicklm)
# scatter plot with regression plane
scatter3D(x, y, z, pch = 18, cex = 2, theta = 5, phi = 20, ticktype = "detailed", xlab = "wing", ylab = "drum", zlab = "price",  surf = list(x = x.pred, y = y.pred, z = z.pred, facets = NA, fit = fitpoints), main = "Chicken Price Analysis")

2.1 First bite

Let me compare the predicted values from the model with the actual values. In other words, look at the residual values.

kfcw$fitval = chicklm$fitted.values
kfcw$residuals = chicklm$residuals
# kfcw
xkabledply(kfcw, title = "Data values and residuals") # better display than the line above

Ah, the residual is most negative with the 40-wing combo at -$1.31. That saves me a can of soda. I should go for that and kill myself.

But if I am wiser, and insist on having a one-person portion:

1. 5-piece wing is at +$0.10 (paying a dime too much!) - Not good.
2. 3-piece drum is at -$0.12, saving a dozen pennies, sounds good.
3. The 3-wing-2-drum deal (at -$0.30) is the best value for me!

2.2 Double dip

Hold on! If we plan to do this regularly, shouldn’t we inspect the percentage savings instead?

kfcw$res_2_p = chicklm$residuals/kfcw$price*100
# kfcw
xkabledply(kfcw, title = "Revised Data values and residuals") # better display than the line above

The new metric here tells me:

1. The 10-drum deal at -4.5% is best percentage-wise.
2. The 3-wing-2-drum deal at -3.76% is next.

To me, putting all these different metrics and criteria together, the 3-wing-2-drum is my winner.

Case closed. Mic dropped.

2.3 Confession

Oh wait… I forgot that I am a pescatarian (vegetarian and fish) for over a decade! I just have these cravings for wings once in a while. What if I add a sin-ful parameter as a penalty term, how much should I order next time to minimize my sin (Root-Mean-Squared-Sinfulness).

To be continued…

But yes, the wings there were really really good. Totally worth the wait.